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COMPUTER SIMULATION OF THE VELOCITY DIFFUSION 


OF COSMIC RAYS 

Thomas B, Kaiser 
Thomas J. Birmingham 
Frank C. Jones 

ABSTRACT 

Monte Carlo simulation experiments have been per- 
formed in order to study the velocity diffusion of charged 
particles in a static turbulent magnetic field. By following 
orbits of particles moving in a large ensemble of random 
magnetic field realizations with suitably chosen statistical 
properties, a pitch-angle diffusion coefficient is derived. 
Results are presented for a variety of particle rigidities 
and rms random field strengths and compared with the pre- 
dictions of standard quasi-linear theory and the non-linear 
partially-averaged field theory described in a companion 
paper. 
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COMPUTER SIMULATION OF THE VELOCITY DIFFUSION 


OF COSMIC RAYS 


I. INTRODUCTION 

It is now widely acknowledged that the velocity diffusion of charged particles 

induced by random fluctuations in a static magnetic field can be fully understood 

1 

only in terms of a non-linear theory. The quasi-line ar analysis of the process, 
which works well in most regions of phase space, breaks down in the case of 
particles whose pitch angles relative to the average magnetic field, 6 , approach 
90°. 

During the past few years, several non-linear theories have been proposed. 
m 2 

The accompanying paper, hereafter referred to as Paper I, gives a detailed 
presentation of one such theory along with a list of references to several others. 
All of them, as well as the quasi-linear treatment, compute coefficients for 
diffusion in p = cos 9 which are mutually consistent for Ip I 1. For p -»■ 0, 
however, they obtain results which vary widely as to their dependence on pitch 
angle, rigidity and random field strength. 

Choice among the competing non-linear theories has been difficult because, 
on the one hand, none of them includes rigorous internal validity criteria and, 
on the other, no experimental observations with which to compare their di- 
vergent predictions are currently available. 

In order to help remedy this situation we have performed computer simu- 
lation experiments for a variety of particle rigidities and rms random field 
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strengths. In this way we have been able to study the diffusion process directly, 
independent of any theoretical assumptions. This has made it possible to clearly 
delineate the region of phase space throughout which quasi-linear theory breaks 
down and to derive the pitch angle diffusion coefficient as a function of n for a 
range of experimental parameters. 

We use the term "experiment" to describe the following procedure. A dis- 
tribution of particles in phase space is allowed to evolve under a particular set 
of initial and boundary conditions in one realization chosen from a previously 
prepared statistical ensemble of stochastic magnetic fields. Exact particle 
orbits in the given field are followed, and information about the time evolution 
of the distribution is stored in the form of various reduced distribution functions 
and their moments. This process is then repeated for each of the several hun- 
dred magnetic field realizations comprising the ensemble, and, finally, the 
results are ensemble averaged. 

Two qualitatively different sets of initial and boundary conditions were used. 
With one set the transient relaxation of an initially anisotropic distribution of 
particles was followed. This design was used primarily as a test of the basic 
simulation method and the computer codes used to implement it. Under the 
other set of conditions a steady state was established by uniform injection and 
absorption of particles. It was from this latter design that we derived our 
principal experimental results. 

The magnetic field model and the method used to generate the ensemble 
are described in Section H. In Section m the time-dependent experiment is 
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discussed. Section IV gives a detailed description of the design and interpre- 
tation of the steady-state experiments, and Section V presents their results for 
a variety of parameter sets. The formulation and performance of the particle 
orbit integration algorithm, a derivation of the appropriate boundary conditions 
for the steady-state experiments, and an analysis of the uncertainty in the ex- 
perimental diffusion coefficient are treated in appendices. 

H. THE MAGNETIC FIELD 

For reasons given in Paper I, each realization of the magnetic field was 
taken to be of the form 

B(r) =^ v 5B(z)+ J eL<B> (1) 

/a/ rv ^ " 

Q 

with 5B(z) a homogeneous Gaussian process of mean zero and two-point 
correlation function 


CgCz'-z) = <6B(z) 5B(z'-)>'/ <5B 2 > 


- exp 


Iz'-zl 
z c _ 


( 2 ) 


which corresponds to the power spectrum 


<P(k) = 



dre-^C B (f), 


(3) 


l+k 2 z 2 * 

V 

The choice of an exponential correlation function ensures that 6B(z) is a 
Markovian process, i. e. , one that is completely specified by its two-point 

4 

probability density. 
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A realization of such a process can be generated on a grid of spacing h by 
the following procedure. 

1. Choose 5B(o) from a Gaussian distribution with mean zero and variance 

<5B 2 >. 

2. Let P [6B(z 7 ) I SB(z) j be the conditional probability density for 5B(z' ), 
given that 6B has the value 6B(z) at z. It is straightforward to show that 
P [5B(z') |SB(z)l is Gaussian with mean SB(z) C b (z' - z) and variance 

( SB 2 } [ 1 - C B (z' - z) 2 ] . Choose SB(h) from a distribution with density 
P [ 5 B(h) ISB(o)]. 

3. Repeat step (2) for succeeding grid points, each time selecting 5B(nh) 
from a distribution with density P [SB(hh) I SB(nh - h)] . Since 5B(z) is 

a Markovian process, the distribution of 5B(nh) depends only on 5B(nh-h) 
and not on 5B(nh - mh) , 2 < m < n. 

This method was used to generate an ensemble of 800 realizations of the 
random component of the magnetic field. Each realization was defined on a 
grid of length 200 z c and spacing h = z c /64 and had <5B 2 ) = 1. Between grid 
points values of the field were defined by linear interpolation. The ensemble 
was stored on magnetic tape and used repeatedly. During the performance of 
an experiment the realizations were read sequentially and the field values scaled 
to give the desired rms random field strength. A segment of a typical realization 
is shown in Figure 1. 

If the empirical correlation function of SB(z) over a finite ensemble of R 
realizations is defined 
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R 


rf fe) = 


<SB 2 >(R-1) 


^ s~ ^ 

« 5B( n )(o)-i- Y 5B( m )(o) 


m = 1 


1 R 

SB( n >(z) -~YL SB {fi) (z) 

R 2= 1 

— *- L J^/ 


( 4 ) 


where the superscripts n, m - , K on the right side index the finite ensemble, then 
it is straightforward to show. that the mean and standard deviation of rj^(z) are 

<rj(z)> = CgCz), (5) 


<[rt(z)-(r^(z»] 2 >^ = 


•t 1 


+ C B (z) 


\r 


R 


Vl 


( 6 ) 


ha the above, the angular brackets indicate averages over the infinite Gibbsian 
ensemble. 

The right side of Equation (4) for the full 800-realization ensemble is shown 
in Figure 2 compared with the theoretical expression. Equation (2). One- 
standard-deviation limits, as computed from Equation (6) are also shown. 


m. THE TIME DEPENDENT EXPERIMENT 

In order to check the validity of our simulation method and the computer 
codes used to implement it, an experiment was designed to monitor the evolution 
of a distribution of particles in a region of phase space where quasi-linear theory 
is valid and the diffusion process, well understood. 

hi this experiment we followed the mean and variance of the space and gyro- 
phase averaged distribution function, f 0 (P, t), as it evolved from the initial 
condition 
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f 0 0,o) = 8Hul-{i 0 ). (7) 

This was accomplished by integrating the trajectory of a single particle in each 
of R realizations of the ensemble of random magnetic fields. Each particle had 
the common starting point fi(t = o) = tx Q and a velocity phase angle chosen ran- 
domly from a uniform distribution on [0, 2 -jt], and was followed for the same 
length of time. The mean and variance were computed as 

, R 

m(t) = <g(t)> R = - £ /i (n) (t) (8) 

R n = l 

1 R 

o 2 (t) = <[M(t)- <M (t) > 1 2 > R = - — - X [/x^(t)-<M( t)>3 2 (9) 

R - 1 |i = ! 


where the superscript (n) indexes realizations over the ensemble. 

The behavior of m(t) and a 2 (t) predicted by quasi-linear theory can be ob- 
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tained from the kinetic equation for f Q 

9f 0 3 f 3f 0 1 

3f c 

6 

where the quasi-linear diffusion coefficient is 


Pgjf o, t) = D Q L (ju, ~) [1 + J(M ? t)] 


( 10 ) 

(ID 


( 12 ) 


DQL^oo) = 


z c <5co 2 ) Iju!(1-m 2 ) 


v 2 n 2 +z 2 ul/v 2 


(13) 


J(/a, t) = e~I ^ I vt/z c [ CO so; 0 t — - — - sin co 0 t] 

ljU.lv 


(14) 
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with 


<5w 2 > = 


q 2 (SB 2 ) 
7 2 m 2 c 2 5 


w o = 


The transient component, J(y, t), is included for completeness. Its effect 
damps out on the correlation time scale, tj?- L = z c (vlft !) -1 , which in the quasi- 
linear regime is comparable to the transit time, z c /v. The contribution of 
J( y, ,t) to m(t) was unobservably small in the time- dependent experiment, and 
while observable, its small contribution to a 2 (t) damped in a few transit times. 
Consequently, in calculating m(t) anda 2 (t) from Equations (10) — (12), J( y, t) 
will be ignored. 


Integration of Equation (10) forward in time yields 


f 0 (ji,t) = 5 (y-y 0 ) + 


3 r 

/ dt' — I 

Jo 0M L 


DgJ (n, -) 


9f 0 (M') 


The mean value of y is then 


m(t) = 


■L 


dy yf Q iy,t) - y 0 + 


f +1 9 rtT 9f 0 (/x,t')~ 

dt' / dy M— Dg£0*,°°) - . (16) 

J- 1 3ju L dy J 


which may be integrated twice by parts using the boundary conditions. Equation 


(11), to give 


m(t) = y Q +J dt' J" dy [f)gg (y, <»)] f 0 (M,t') 


Furthermore, in the region of y -space covered by this experiment, the variation 
of Dgy(ju, 00 ) is very nearly linear 
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DJJH^OO) S Dg L +Dj QL (M-Mo), 


(18) 


Upon using this form plus the normalization 

dpf o 0i,t') = 1 

we obtain 

m(t) = M 0 +Dp L t (19) 

The variance of f 0 , 

o 2 (t) = f dp [j£-m(t)l 2 f 0 (/x, t), (20) 

•'-1 

can be derived in a similar fashion by computing the second moment of Equation 
(15). The result is 

a 2 (t) = 2 D^ l t + (Dp L t) 2 (21) 

Data for the ease z c co 0 /v = 1, z^Sco 2 ) 1 ^ 2 / v = 10 -3 / 2 , R = 200, p 0 = (3/4) 1 / 2 
are given in Figures 3 and 4 with Equations (19) and (21) plotted for comparison. 
Note especially that the transient component of has negligible effect on m(t), 
and its effect on a 2 (t), which is responsible for the initial depression below the 
value given by Equation (21), damps completely in a few transit times. The 
excellent agreement between the behavior predicted by quasi-linear theory and 
that actually observed in the computer simulation strongly indicates that the 
simulation code accurately depicts the physical phenomena. 
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XV. THE STEADY STATE EXPERIMENT: DESIGN AND INTERPRETATION 


a. Experimental Design 

In order to study the diffusion process in regions of phase space where 
quasi-linear theory breaks down, experiments of the type described in Section 
HI were found to be inadequate. Transient effects damp so slowly that by the 
time they are negligible the initial distribution function has spread to the point 
where Equation (18) is badly violated. In other words, the correlation time 
becomes comparable to the relaxation time of the initial distribution. An illus- 
tration of this behavior is given in Figure 5, in which a 2 (t) as computed from 
Equations ( 8 ) and ( 9 ) is plotted for the ease z c cj 0 /v - 1 , z c < 5 uj 2 ) l ! 2 lv = 10 " 3 ^, 

R - 200 and m 0 = 0. The initial transient oscillations are seen to persist for 
more than 50 transit times and clearly are not negligible. By the time they have 
apparently damped away, the distribution has relaxed to a point where its evolu- 
tion is no longer characterized by the value of the diffusion coefficient at g = p 0 . 

An obvious way to circumvent this difficulty is to design an experiment in 
which the relaxation time is very long, in fact infinite, i. e. , one in which the 
particle distribution has evolved to a steady state. 

Of course, the only steady-state solution of a diffusion equation like Equation 
(10) with boundary conditions Equation (11) is the uninteresting isotropic distribu- 
tion, f Q = constant. This situation would be altered, however, if a source of 
particles and absorbing boundaries were introduced. Then, when the influx of 
particles from the source was balanced by the escape of particles through the 
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boundaries, a steady state would be attained. Physically, such a situation would 
be similar to the injection of particles into a magnetic mirror field and their 
subsequent escape by diffusion into the loss cones. 

b. Computer Implementation 

The steady state experiment was executed in the following way. A total of 
N particles with independent and uniformly distributed velocity phase angles 
were launched at t = 0 at the point :r = r s with ju = p s in each of R realizations 
of the magnetic field. Typically N was 10-20 and R was 800. Absorbing bound- 
aries were located at p = < fi s and ii = p R > /i s . 

The trajectory of each particle was followed until it left the region 
H L ix R , at which time the particle was annihilated. Snapshots of the 
evolving particle distribution function in each realization were taken at regular 
time intervals At for a total time t F , long enough for all particles to be absorbed. 
The R collections, of K = t F /At snapshots each, thus generated were then aver- 
aged over the R realizations of the magnetic field to give an ensemble averaged 
collection of snapshots, <G) r (ju, t k I n s , 0), t k = kAt, k = 0, 1, . . . , K. 

We show now that 

K 

<f> E - X) <G> R (e,t k > s ,0) 

k - l 

= -jf dr(G> R (ju,rip S3 0) (22) 

satisfies approximately the steady-state diffusion equation 
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Djum^’ 00 ) 


d<f>. 


q- 5 . 

J At 


S(m-m s ) 


(23) 


dfl l ^ djLl 

and equivalently that <f> R is the steady-state one particle distribution function 
averaged over the finite ensemble. 

To justify this assertion and approach, consider <G> r (m> t f m s * t ' ), the en- 
semble averaged Green's Function corresponding to an injection source of 
strength N located at '{jl = at time t\ <G> R satisfies the equation 


a<G>, 


R 


at 


a 

a 11 


(Mj t - 1') 


a<G> fi 


3ja J 


+ N5(ju-M s )8(t-t') 


(24) 


the initial condition <G> R (ju, t i ,u g , t') = 0 for t < t' and the boundary conditions, 
derived in Appendix B, 


, , a<G> R 

D (M CO) * ± <0-)<G> r 

MM a R 


= 0 


(25) 


M=M R ,M L 


where the upper and lower signs apply at the absorbing boundaries at m r and m l 
respectively and (a) is a m - dependent quantity given in the Appendix. The 
"turn-on" time of the diffusion coefficient is taken to be the injection time 
t' , since fluctuations, 5G , x from Realization to realization in the ensemble 
commence growth from zero at this time (cf. Paper I). 

Note that since time enters Equation (24) only in the form t - t 1 , '(G) R is 
time translationally invariant 


<G> r (m, t t fi s , t') = <G) r (m, t -t'U s , 0) 


11 



Consider then the quantity 


f dr<G> R (m, tl/u s ,r) = f dr <G> R (ft r I m s , 0) 

*'o *'o-’ 

•7 

where, as in the usual Green’s Function prescription, the integration extends 
to t + = t + e. The relationship between <f) R compiled in the computer experi- 
ment and the solution to Equation (24) is quite obvious 

1 f X+ 

< f >R = t S2o <F(t)) * = t^ At J Q dT<G> R (/i ’ tIp S’ r ) 

To verify Equation (23), note that differentiation of <F(t)> R with respect to 
time and use of Equation (24) and its initial condition result in 


0 < f> r i r t+ a<G> R (p,ti Ms ,r) 

= — / dr 

3t At Jr, 3t 


= ~ f d r^~ [d ( ft r) 

At Jo L W 


9<G> r (ft rim, 0)1 N 


+ — 5(p-n s ) 
At S 


So long as D^(ft t) is asymptotically positive, i.e. , that a time t A exists 
such that (ju, t) > 0 if t > t A , then an H-theorem for <G> R is easily derived 
from Equations (24) and (25) by multiplying Equation (24) by <G> R and integrating 
over p. Consequently <G> R ( ft r ija s > 0) -* 0 as ’r -»■<» so that <F> R approaches a 
steady state. Thus, asymptotically in time. Equation (26) reduces to 


[f"> 


3<G) r (/i, tI/Jc, 0) 

,(M,r) £ 5 = - NS(m-« s ) 

dii 


The diffusion coefficient D^(/x,r) saturates, i.e., approaches the asymp- 
totic value D (n, ”), in a time of the order of the correlation time r c s z c /(ijulv) 
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of the fluctuating magnetic field (cf. Paper I). There is also a finite transport 

time of particles from the source point m s to the point Roughly speaking 
% I p - i 

r t _ 2 — - , so that <G> R ss O for r < r tr (ju) and this region of small r 

1 

does not contribute to the integration in Equation (27). 

Provided that 

t c (m) < r tr (/x) ( 28 ) 

we may then write 

(29) 

= M IW(M '“ ) ^ <f>R ] = ‘S 5(f -^ 

Note that in the immediate vicinity of the source Equation (28) is violated and 
hence Equation (29) is invalid. This violation is manifest, as we shall see, in 
the computer experiment. In our analysis of the computer results we shall 
quantify the region of validity in ( m - M s ) of the approximate form Equation (29). 
Note from Equation (27) that the particle fluxes 



dr D p/M T ) 


8<G) r (ju 3 t|p s ,0) 

3p 


to the left and right of the source are constant, independent of ju . Thus in the 


region where Equation (29) is valid 


- J T - — ='- J R 
L At R 


H L < [i <li s 

[x s <[i<n R 


where J L and J R are the constant values of the left and right fluxes. 
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j /-C\ 

In the computer experiment . — d R - , J L , and J R are directly and separately 
measured so that D (pi, °°) can be determined readily from Equation (30) in- 

fAJJ, 

dependent of any theoretical assumptions save that Equations (24), (25), and 
(28) are valid. 

c. Results of a Typical Experiment 

In practice, <f> R (pi) was obtained in the form of a histogram of bin width. 
Apt , typically = 10 -2 . The result for the case z c co 0 /v = 1, <5 co 2 > ^ /m 0 = 0. 1, 
pt L =-0.65, pi g = 0.405, pi R =0.65, Api = 0. 01, ( v/z c )At = 0. 03125, N = 8, 

R = 800, which was typical of the parameter regime throughout which simulations 
were performed, is given in Figure 6. 

Before proceeding to a more detailed analysis of experimental results , we 
briefly discuss various gross features of the steady state distribution evident in 
Figure 6. Perhaps the most obvious characteristic of the distribution is the - 
presence of particles on the opposite side of pt = 0 from the source. The appar- 
ent ease with which particles scatter through 6 = 90° is in direct contradiction 
to the prediction of quasi-linear theory, which yields a vanishing diffusion co- 
efficient at pi = 0. 

Also evident is a flattening of the distribution function just to the left of the 
source. Although the diffusion coefficient and, therefore, according to Equation 
(30), d <f> R /dpi, are expected to be even functions of pi, such is clearly not the 
case. We ascribe the lack of symmetry to the breakdown of the adiabatic ap- 
proximation near the source and estimate the extent of the effect in Section IV-e. 
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Finally, we call attention to the fluctuations in the number of particles 
occupying each histogram bin, the standard deviations of which comprised one 
of the diagnostics available in the experiments. It was observed that the magni- 
tudes of these fluctuations, relative to the mean bin occupation numbers, were 
of order = 0. 04, as would be expected on the basis of simple statistical 
arguments. This scaling gives an accurate measure of the' cost entailed in im- 
proving the statistical accuracy of the experiments. 


d. Determination of (g, °°) 

The diffusion coefficient, (p, 00 ), was obtained from Equation (30), 
using the observed particle fluxes and the measured slope of the histogram 
representation of <f> R (ja). 

The particle fluxes were determined by counting the number of particles 
that escaped through each boundary. Thus, 


J 


L 



(31) 


J 


R 



(32) 


where N L and N R are the' ensemble averaged number of particles that escaped 
at (i = n L and M = M R , respectively. 

The slope at a point, p , was obtained from a linear least squares fit to the 
histogram points symmetrically bracketing g. Typically, ten point’s were used, 
the actual number being determined by requiring the - statistical uncertainty in 
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to be -10%. A detailed analysis of the -uncertainty in D^is given in Appen- 
dix C. 

The diffusion coefficient computed in this way from the data in Figure 6 is 
given in Figure 7. Also shown are the quasi-linear expression. Equation (13), 
and the diffusion coefficient derived in the partially averaged field theory de- 
scribed in Paper I, 

Several details of the p-dependenee of the experimental diffusion coefficient 
apparent in Figure 7 proved to be generally characteristic of D JUJU ( p, °° ) through- 
out the parameter regime studied, and we call attention to them here. 

First of all, to within statistical uncertainty, D MjU (p, “) is an even function 
of p. This fundamental property, predicted by all theories, is of interest chiefly 
as a check on the simulation scheme. 

Furthermore, there is a significant peak in D MM (p, °°) at p = 0. The en- 
hanced scattering represented by this peak is predicted by the partially averaged 
2 

field theory, which attributes the effect to particle mirroring. 

Finally, the pitch angle regime throughout which quasi-linear theory fails 
is roughly characterized by I p I % < 5 B 2 ) Vz /< B >. A theoretical explanation for 
this result is given in Paper I. 

e. Experimental Justification of Approximations 

There were several approximations necessary to the derivation of the funda- 
mental equations in terms of which we have interpreted the results of our computer 
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experiments, viz,, Equations (24), (25), and (28). The experiments" themselves 
provided information that made it possible to check the validity of these approxi- 
mations a posteriori. We turn now to a consideration of three of them. 


1. Gyrotropy 

The derivation of Equation (23) assumes that the particle distribution is very 
nearly gyrotropic. There were in fact’ small departures from gyrotropy in the 
measured <f> R . To check the validity of the gyrotropy approximation, the 
phase angle distribution function, 


was Fourier analyzed 





h(0) = ^ h n e in< £ , (33) 

n = _ oo 

for each experiment. A histogram of h ((f)) for the of Figure 6 is given in 
Figure 8. The departure of h from h 0 (dashed line) is evidently concentrated in 
the n = ±1 modes and of small magnitude. The actual amplitude of the modulation 
in this case was measured to be 


I h - h Q I 


max 


h 


8.6 x 10" 3 , 


which is in good agreement with the estimate derived in Paper I, viz. , 


I h - In, I (8u> 2 ) 

^ <~w 

h o OJ l 


= 10' 2 
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2. Boundary Conditions 

In deriving the boundary conditions, Equation (25), assumptions were made 
concerning the symmetry and statistical correlations of the probability distribu- 
tions of <p and 6 to. To verify that these boundary conditions were actually satis- 
fied in the experiments, we considered the equivalent Equations (B9), (BIO), 
since each quantity appearing in the latter equations was directly determinable 
from the experimental data. 

Thus, with J and (f > R determined in the manner described in Section IV-c, 
and rv-d, and (a) obtained from Equation (BIO) and the relevant experimental 
parameters, the comparison for the case of Figure 6 is summarized in Table I. 
The close agreement between J <f) R and (a) at each boundary indicates that the 
boundary conditions. Equation (25), are indeed satisfied in the experiment. 

3. Adiabatic Approximation 

As already mentioned, a diffusive description of the behavior of <f ) R in the 
vicinity of the source appears to be invalid. Consequently, Equation (30) cannot 
be used there to determine the diffusion coefficient. In Section III-B, we dis- 
cussed this difficulty from the Green's Function point of view. We here discuss 
the same breakdown from an alternative but equivalent. approach. 

It is shown in Paper I that within the context of the partially-averaged field 
theory, the general steady state equation satisfied by f 0 , the gyro-phase averaged 
distribution function averaged over the infinite ensemble is 
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(34) 


l{fs <>■•>,), 

where <Af 0 is the asymptotic (in time) change in f 0 apparent to an observer 
moving from the point pt at t = 0 along an orbit in the partially- averaged field 
(SB)^ , and < denotes a final ensemble average over the distribution of 
( 83 )^. 

If f 0 is expanded in a Taylor series about (x , this change can be represented 

9f 

<Af 0 V = + 0«Afi> 2 f ) (35) 

so that Equation (34) becomes 



= -(^) 8(C-^) 


which is of the same form as Equation (23) if higher order terms are neglected. 
Within the context of the partially-averaged field theory, use of Equation (36) 
neglecting 0«Ap>^) terms to approximate Equation (34) is an adiabatic approxi- 
mation similar to that made in standard quasi-linear theory. 

Since the neglect of the higher order terms in expansion (36) is clearly 
invalid near the source, where 3f 0 / 3/u is discontinuous, a diffusion equation 
like Equation (36) should not be expected to apply there. 

Integration of Equations (34) and (36) from pt = pi L to an interior point re- 
sults in two different expressions for the flux: 
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In regions where a diffusion description is valid these expressions should agree 
with each other and with the right side of Equation (30). J non ad and J ad , will 
disagree where a diffusion picture, and thus Equations (36) and .(30), are in- 
appropriate. For a given experiment these regions can be mapped by computing 
J non ad and J ad using the experimental distribution function <f) R for f Q , and 
comparing. Such a comparison is shown in Figure 9 for the case of Figure 6. 

It is evident from the figure that Equation (30) cannot be used to determine 

00 ) in the region 0. 2 < p < 0.65. Also shown in Figure 9 is the flux 
observed in the experiment. The discrepancy between the computed flux and 
the observed flux is due to the fact that the partially-averaged field theory 
slightly overestimates the scattering at i± - 0, as is clear from Figure 7. 


V. RESULTS 

Steady state experiments like the one described in detail in Section IV were 
performed for a variety of values of the dimensionless parameters, 

Z C W 0 

6 = , 

V 


which measures particle rigidity, and 

(6m 2 ) V. Vl 

V = > 

w 0 
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which characterizes the strength of the random magnetic field. For the case of 
Figure 6 , for example, these values were e = 1, rj = 0.1. 

Results of the experiments are given in Figures 10 a-h, where the dimen- 
sionless diffusion coefficient, (z c /v)D^ (p, 00 ), is plotted as a function of M 
for each of the ( e, rj ) combinations tried. Also shown for comparison are the 
predictions of standard quasi-linear theory and partially-averaged field theory, 
represented by dashed and solid curves, respectively. The experimental pa- 
rameters used in each of the experiments are given in Table II. 

As mentioned in Section IV-E, a velocity phase angle distribution function 
was compiled for each experiment to check the assumption of gyrotropy. In 
general, departures from cylindrical symmetry were concentrated in the n = ±1 
modes and had amplitudes in accordance with the theoretical estimate of Paper I, 


h-h 0 l 


= V 2 . 


(39) 


A plot of max {lh - h Q l}/h 0 as a function of 77 is given in Figure 11, with the 
estimate (39) superimposed for comparison. In no case was the departure from 
gyrotropy greater than 5%. 

Since the disagreement among the various theoretically computed diffusion 
coefficients is greatest at 11 = 0 , it is useful to investigate the manner in which 
CO ) scales with e and 77 . The scaling observed in the experiments over 
roughly one decade of variation in each parameter is given in Figures 12a, b. 
From the figures it is clear that D^(0, 00 ) has a power law dependence on each 
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parameter, with the exponent an increasing function of the other parameter, i,e. , 


D /xju (0,«)ai?m( e ) e<lOrt . 

Least squares fits to the data give m(0. 25) = 2. 5, m(l. 0) = 2. 7, m(4. 0) 
and q(0. 05) = 0.18, q(0.10) = 0.33, q{0.30) = 0.62. 


(40) 
= 3.1, 
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Table I 


Data Needed to Check Satisfaction of Boundary Conditions for the Case of Figure 6 


M 

(z c /v)l J l 

<f> R 

(z c /v)l Jl /<f> R 

(z c /v)<a> 

-.65 

+36.4 ±1.1 

854. 3 ± 62.8 

.043 ±.003 

.039 

+.65 

219.6 ±1.1 

5319. 1 ±137.1 • 

. 041 ±. 001 

.039 



Table H 


Experimental Parameters Used in the Steady State Experiments of Figures 6 and 10 


Exp. 

6 

n 

R 

N 


% 

^ R 

Aju 

vAt/z c 

1 

0.25 

0.05 

400 

10 

-.65 

.405 

.65 

.010 

.25 

2 

0.25 

0.10 

400 

32 

-. 65 

.405 

.65 

.010 

.25 

3 

0.25 

0.30 

800 

20 

-.65 

.405 

.65 

.010 

.25 

4 

1.00 

0.05 

100 

80 

-.65 

.405 

.65 

.010 

.0625 

5 

1.00 

0.10 

800 

8 

-.65 

.405 

.65 

.010 

. 03125 

6 

1.00 

0.30 

800 

8 

-. 70 

.605 

.70 

.010 

. 03125 

7 

4.00 

0.05 

800 

8 

-.50 

.3025 

.50 

.005 

.015625 

8 

4.00 

0.10 

800 

10 

-.65 

.405 

.65 

.010 

.015625 

9 

4.00 

0.30 

800 

16 

-.35 

.4025 

.55 

.005 

.00390625 
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8B(Z) (ARBITRARY UNITS) 



Figure 1. Segment of a. typical realization of the random component of the magnetic field. 


CORRELATION FUNCTION 


1.10 



Figure 2. Empirical correlation function, of 5B(z) over the full ensemble of 800 realizations, Equation 
(4). The theoretical expression. Equation (2), is also plotted. 


25.00 




Vt/z c 

Figure 3. <J“ (t)> in a time dependent experiment with z c wjy = 1, z c = 

10" 3/2 , R = 200, m 0 =i(3/4) 1 ' 4 . Equation (19) is plotted for comparison. 
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7 

Figure 4. ([mW - <m(*) >I 2 > in the experiment of Figure 3, with Equation (21) 
plotted for comparison.- 
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0 1 2 

log (vt/z c ) 

Figure 5. Same as Figure 4 except m 0 = 0. 


30 



NUMBER OF PARTICLES 



Figure 6. 


Histogram of f 0 ( n , «°) in a steady state experiment with z c co Q /v = 1, <5 co 2 } ”/ co Q = .1, 
M l = - -65, = -405, /u R = .65, A ^ = .01, v At/z c = 2“ 5 , N = 8, R = 800. Error bars' are 

one-standard-deviation statistical uncertainties. 



(Zc/V)*D(MU,MLj) 









4 > (DEG.) 

Figure 8. Histogram of the velocity phase angle distribution function, h ( defined by Equation 

(33) for the experiment of Figure 6. The gyrotropic component, h 0 , is indicated by the 
dashed line. 



Figure 9. 


Flux in the experiment of Figure 6. Sho 
the observed J, and J„ . 



areJ nonad Equation (38), J ad Equation (39), and. 






<Zc/V)*D(MU/MLl) 



Figure 10a. Dimensionless diffusion coefficient, ( 36 /v) D ^ derived from steady state experi- 
ments whose parameters are given in Table II. Also shown are the predictions of partially- 
averaged-field theory (solid line) and standard quasi-linear theory (dashed line). 




(Zc/V)*D(MU/MU) 


3.00T Xl0f-3 


EXPERIMENT 3 



Figure 10c. 


0.75+ 


/V)*t>(MU/MU 







(Zc/V)*D(MLJ/MU) 





(Z</V)*D(MU/MU) 



MU = C05CTHETR) 


Figure 10h. 


0.7S- 



-2 -1 0 


log r] 

figure 11. Plot of max j lh - h 0 ! [/h G , which measures departure from gyrotropy, 
for each of the experiments of Table IT. The theoretical estimate. 
Equation (40), is superimposed for comparison. 
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log [(z c /v) D^(O,00)] 


-1 



-1.6 -1.4 -1.2 -1.0 -0.8 -0.6 . -0.4 -0.2 0 

log 17 


Figure 12a. Dimensionless diffusion coefficient at fx = 0, (z c /v) D derived 

from the experiments of Table II. a) (0, <») as a function of 17 for 
constant values of e . Solid lines are least squares fits. 
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log[(z c /v) D /i/i (0,oo)] 





.8 .6 .4 .2 0 .2 .4 .6 .8 

log £ 


Figure 12b. Dimensionless diffusion coefficient atj i = 0, (z /v) D (0, “), derived 

Mr 

from the experiments of Table II. b) D^ M (0,“) as a function of e' 
for constant values of r? . Solid lines are least squares fi-fcj. 
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APPENDIX A 


THE PARTICLE PUSHER 


I. THE ALGORITHM 

The numerical scheme for integrating the particle equations of motion ex- 
ploits the fact that for time increments, t' - t, satisfying Igj (x) (t' - t)| <1, 
where to(x) is the local gyrofrequency, the effect of a static magnetic field on a 
particle's velocity is simply to rotate it about the total field. Thus, 

v(t') = T(t, t') v(t), (A.1) 

where the transformation matrix, T, is a rotation. 

rv 

Propagation in configuration space is accomplished by means of the second- 

O 

order implicit algorithm 

x(t') = x(t) + [v(t) + v(t')j . - (A.2) 

rw aj aj j 


The remainder of this section is devoted to deriving the elements of the 

matrix T. Consider two Cartesian coordinate systems, S and S', whose z-axes 
fa 

are aligned respectively along the average and total magnetic fields. If a unit 
vector along the total field is denoted b, and the coordinate unit vectors, 
el, then S' is defined by 


e 2 ' - (t 3 x S)/le 3 x Si 
e/ = (e^ X S) X S/l^ X Si 


A-l 



In S' the velocity vector evolves according to 



v'(t') = gv'(t), 
cos co(t' - 1) sin co(t' - 1) 0 

-sin toO' - 1) cos c«j(t' - 1) - 0 


0 0 1 


(A.3) 


(A.4) 


The gyrofrequency, of course, depends on x. It will be seen below that if the 
algorithm is to be reversible in time and maintain second order accuracy it is 
necessary that the position at which co is evaluated be the intermediate point. 


XjCt, t') = x(t) + - v(t) (t' - 1) 


(A. 5) 


The transformation connecting S and S' is also a rotation: 

v' = R v (A. 6) 

In terms of the magnetic field components, R is 


R = 

where B = IB I and B^ = IB - (B * e 3 ) e 3 1. 

If Equation (A. 6) is used in Equation (A. 3) and the result compared with 
Equation {A. 1) it is clear that 


B 1 B 3 

BB, 


b 2 b 3 

BB, 


Bjl 

B 


B, 


B, 


0 


(A.7) 


B, 


Bo 


Bo 


B 


B 


B 


T = R _1 R 




(A.8) 
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In terms of a vector, co, 


co = a) b 


and the functions 


*1 


sin co(t' - t) 1 - gos cb(t' - 1) 

] ■> %2 ~~ ^ 2 
co(t'-t) oo 2 (t'-t ) 2 


T has the explicit form 


T n = + 

T 12 = ?i « 3 (t'-t) +'-| | 2 co lt o 2 (t'-t) 2 
T 1 3 = C0 2 (t'-t) +| g 2 CO^t'-t) 2 

T 21 = ~h J g 2 Wl « 2 (t'-t ) 2 

T 22 = 1-^ 2 K + C 0 2 )(t'-t ) 2 

^23 — ^1 W l ^ + ~2 W 2 W 3 ^ 

T 31 = h "2^' -t) + { ^2 C0 1 C0 3 (t' - t ) 2 
T 32 = ^ g 2 C0 2 C0 3 (t'-t ) 2 

T 33 = l-i? 2 (co 2 + co 2 )(t'-t ) 2 


(A.9) 


(A. 10) 


(A. 11) 
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Although apparently complicated, the right sides of Equation (A. 11) are well 
suited for numerical computation because the functions £ 1? £ 2 can be expanded 
in rapidly convergent power series in the small quantity co(t' - t). Consequently, 
the transformation involves only the operations of multiplication and addition; 
no quotients, square roots or other time-consuming operations need be 
performed. 

The particle pushing scheme is completely specified by Equations (A. 1), 

(A. 2), (A. 4), {A. 5), (A. 7) and (A. 8). Various properties of the algorithm are 
described in the following section. 

II. PROPERTIES OF THE ALGORITHM 

a. Time Reversibility 

Since the exact equations of motion of a charged particle in an electromag- 
netic field are reversible in time, it is desirable that the numerical integration 
scheme also possess this properly. This is especially true in the present ap- 
plication, since the effect of irreversibility is a spurious diffusion of particle 
orbits. 

To see that the algorithm is time-reversible, consider the effect of applying 
it twice in succession. A push from t to t' is followed by one from t' to t" = t. 
The phase space position undergoes the transformations ^ v) -* (x' } V) 

(x" , v" ), Time reversibility requires that (x" , v" ) = (x, v). 
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According to Equation (A. 1), v" is related to v by 





= T[xj(t' T [xj (t, t');t, t']y. 


where the dependence of T on Xj has been indicated explicitly, 
tween x" and x is obtained from Equation (A. 2): 


. (A.12>- 
A relation be- , 


x" * x' + ~ (v' + v")(t-t') 

= X + | (V + v'Xt' - 1) + J (V + v")(t - 1') (A. 1 3) 

= £ + \ (v - v")(t' - 1) 


From Equations (A. 12) and (A. 13) it is clear that the scheme is time reversible 
if 


T[x I (t' f t);t',t] T [Xj(t, t');t,t'] = / I 


(A. 14) 


Use of Equations (A. 2) and (A. 5) shows thatXj (t', t) = x T (t, t r ), i.e. , that 
the forward and backward transformations both evaluate the magnetic field at 
the same point. This fact along with Equation (A. 8) reduces Equation (A. 14) 
to 


R" 1 (x T ) SI (t - 1') Cl (t' -t) R (x T ) = I 

« - 1 « « » ~ I fa 

which use of the definition (A. 4) shows to he identically true. Therefore the 
algorithm is, in fact, exactly reversible in time. 


b. Conservation of Energy 

Another property of the exact equations of motion of a particle in a static 
magnetic field is the conservation of kinetic energy. That energy is conserved 
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by the numerical scheme follows immediately from the fact that the kinetic 


energy is proportional to 1 v| 2 , which is invariant under the orthogonal trans- 
formation, T. 


c. Accuracy 

An important consideration in any numerical approximation is its accuracy. 
In this section it will be shown that the accuracy of the particle pushing algorithm 
is of second order in the time step. 

Let x = x(t), x f = x(t'), v = v(t), v' = v(t') and At — t’ - t. Then 

dv 1 d 2 v 

v' = v '+ — At + At 2 + 0(At 3 ) (A. 15) 

~ ~ dt 2 dt 2 

In terms of the matrix 


A(_x) = 


0 a>*(x) -co 9 (x) 

-CO^(x) 0 CO. (x) 

J f-\J *■ 

C 0 2 (x) -OJj (x) 0 


the particle acceleration is 


dv 

-f = Av 
dt 


With this, Equation (A. 15) can be put in the form 


(A. 16) 


(A. 17) 



rv 


v + 


(a+M — ) vAt +>- A 2 vAt 2 + 0(At 3 ) 
W dt 2 / ~ 2 ps 


(A. 18) 
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Now the value of A at the intermediate point, x, = x + — vAt, is 

1 9 'n-I «« Q 


. 3A At dA At . 

A(x t ) - A(x)'+ -a.v — +0(At 2 ) = A + — — + 0(At 2 ). 
3x^2 » dt 2 


Furthermore, 


A(x,) 2 " A(x) 2 + 0(At). 

Therefore, Equation (A, 18) can he written 

t = X + A &i)X At + ^ 4^i> 2 X At2 + °( At3 > (A. 19) 

Since, by the definitions (A. 10), 

M*i) = S 2 (&> = 1 + °( At2 )> 

it is valid to write Equation (A. 19) as 

I = ft+Si^A^^At+l^tx^Afx^At 2 ] v + 0(At 3 ) (A.20) 

Reference to Equation (A, 11) shows that the jnatrix in square brackets in 
Equation (A.20) is precisely T. Therefore, the truncation error of the velocity 

Jv 

integration algorithm is third order in the time step. 

It is easily seen that the same is true in position space. The Taylor series 
for x' is 

dx 1 d 2 x „ 

x' = X + — At + At 2 + 0(At 3 ) 

~ ~ dt 2 dt 2 


At dv ‘ At . 

- x + v — + (v+— At) — + 0(At 3 ) 
~ ~ 2 ^ dt 2 


= 2. + <> + l') -y + 0( ^ At3 ) 

which is identical to Equation (A, 2) through second order terms. 
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d. Computational Speed 


In practice, the functions ^ and £ 2 are expanded in power series in co At 
up to and including second order terms. This maintains second order accuracy 
in x and v and time -reversibility and energy conservation to 0[(coAt) 5 ] . 

In the special case of a slab model magnetic field, where 

co(x) = % 6oj(x • %) + ^ 3 co 0 , 

the computation time required by an IBM 360/91 computer to push a particle’s 
phase space position (x, v) through one integration step is 25/isec. This in- 
cludes a table look-up with linear intepolation to evaluate co^x) and accumulation 
of the velocity distribution function. A typical simulation experiment involving 
10 4 particles, followed for an average of 2000 time steps each, thus requires 
about 8 minutes of cpu time. 
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APPENDIX B 


BOUNDARY CONDITIONS FOR THE STEADY STATE EXPERIMENT 


Let a = dp/ dt and T ( p , a ) be the joint probability density for a particle to 
be found in the incremental range {(p, a), (p + d p, a + da)} . Then the flux 
of particles in p-space is 


J(m) = 



da a SI> (p, a) 


(B.l) 


In order to relate T to f c (p), the ensemble averaged distribution function, 
we also define the conditional probability density for a, given p, \f/(a Ip). Thus 


^(ju, a) = f Q (p) \Kalp), 

so' that 


(B.2) 


J(M) = f 0 (M) 



da a ij/(alp) = f Q (p) <a(p)> 


(B.3) 


When p is an absorbing boundary p = p fi , the conditional probability, 

I Mg )* can be estimated as -follows. In a magnetic field of the form (1), 
Newton's equation for. a is 


a = — — - (1 - ft 2 8 co sirup (B.4) 

dt 

where 5co = qSB( 7 mc) -1 and <p is the gyrophase angle (- it < <p < ir)» For def- 
initeness, consider a boundary to the right of the source, p B = p R > p s . Then, 
since the boundary is perfectly absorbing, there can be no particles with a<0, 
in which ease Equation (B. 4) implies that 6co and <p must have opposite algebraic 
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2 

signs. Furthermore, the distribution of Sco (regardless of 0) is Gaussian, 

and, hence, symmetric about 5co = 0, and the distribution of <f> (regardless of 
2 

6 co ) is uniform, and therefore symmetric about 0 = 0. Assuming that this 
symmetry holds for the conditional distribution of Sco (given 0) and 0 (given 
5 u> ) , i. e. , that for every particle with ( 5 co, 0 ) there is one with ( - Soo , - 0 ) , 
we can equivalently take 

«“ (1 -#*£)*** (B.5) 

where 


X = I Sul, o< y . 


a = I sin 0 1 , 0< <j < 1 . 

Then, with the further assumption that x and a are statistically independent, we 
deduce that their joint probability density is 


<l(X, a) - 2(27r< 6 co 2 exp ( ) — ( 1 - a 2 ) - * 4 

V 2 (6o)2>/ 7T 


(B.6) 


The variable a is constant along hyperbolic contours in the ( x » o ) plane 
according to Equation (B. 5). The total probability that a is less than a given 


value, , is, therefore 


P {*<%}=[ da f 

* / o • / o 




dx q(X, a) 


(B.7) 


The probability density for a is 

3p{«<«o} 


0(«Im r ) = 


30in 


a 0 = a 


(B.8) 


1 r a 

dp - 1 ■ ■ - ;s . 17 

a j_a(l - lit, T J 
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with the function, q, given by Equation (B. 6). Although the integration over a 
can be done analytically,, reference to Equation (B. 3)' shows that only the first ’ 
moment of \p(a\ fx R ) is needed. Thus, when Equation (B. 8)- with Equation (B. 6) 
is substituted in Equation (B. 3) and the double integral evaluated, the result is 

J (%) = «M R )> f 0 (M R ) (B.9) 

fa(n K )> = (|) 3,2 (du 2 } 1 ' 2 (1 -h 2 r ) u1 . (B;10) 

In a system described by a diffusion equation like Equation, (30), conserva- 
tion of particles requires that the particle flux be given by 

i<m) = - •w 5- (bji ° 

Combination of Equations (B, 9) and (B, 11) then gives the boundary condition on 
£ o: 

= 0, M r >M s . (B.12) 

M "* 

Equation (B. 12) was derived for the ease pi R > m s . A similar derivation for 
Ml < M s results in the boundary condition 

= 0,p L <p s . (B.13) 

m = m l 

In the presence of absorbing boundaries the Green's Function <G) R also satisfies 
conditions (B.12) and (B.13). 
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These boundary conditions have a simple physical interpretation. If the 
diffusion coefficient is used to define a "mean free path" for scattering in- /*- 
space as 



= iaUi ))- 1 , 

then Equations (B.12), (B. 13) can be written 


f 3f o 1 "1 

° ± — f 

X °J 


(B.14) 


(B.15) 


Equation (B.15) says that f c vanishes at a point one mean free path outside the 
absorbing boundary. In this form the boundary conditions closely resemble 

7 

those given by Morse and Feshbach in an analysis of spatial diffusion near an 
absorbing boundary. 
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•APPENDIX C 


ERROR ANALYSIS 


The uncertainty in the value of the diffusion coefficient derived from .Equation 
(33) is a composite of .the statistical uncertainties in J and d<f) R /djti. The en- 
semble averaged flux is 


«> - { t 


(C.1) 


£ = 1 


where is the flux observed in the £-th realization. The uncertainty in <J> 


is taken to be 


t [J*- «>]’}“• 


(C.2) 


Similarly, if is the number of particles in the i-th histogram bin (width Aju, 
centered at jUj) in the £-th realization, the ensemble averaged distribution 
function is 




(C.3) 


£ = l 


with uncertainty 


iTT,?, 


% 


(C.4) 


Typically, a 3 /(J> and /<f> R (*ij) were a few percent or less. 

The gradient of <f > R at ju, d<f ) R /dju, was taken to be the coefficient of the 
linear term in a weighted linear least-squares fit to the ten or so points. 
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(f ) R (juj), symmetrically bracketing p. The uncertainty in d<f> R /dp is then 

9 

given by the standard formula 


E4 


a d<f > R /dp , 2 p „ -i 2 

e-^eMe§ 


'fj i °\ LV ff f.J 


The resultant uncertainty in ( p, co )> then, is 




a 


d<f> R /d M 


Dpp L <J)2 (d<f> R /dM) 2 J 


% 


The number of points actually used in the fit was adjusted to give a D /D MjU 
provided that this did not require an unacceptably large spread in p. 


(c: 5) 


(C.6) 

. 10 , 
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